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ABSTRACT 

We have studied non-axisymmetric standing accretion shock instability, or 
SASI, by 3D hydrodynamical simulations. This is an extention of our previous 
study on axisymmetric SASI. We have prepared a spherically symmetric and 
steady accretion flow through a standing shock wave onto a proto-neutron star, 
taking into account a realistic equation of state and neutrino heating and cool- 
ing. This unperturbed model is supposed to represent approximately the typical 
post-bounce phase of core-collapse supernovae. We then have added a small 
perturbation (~1%) to the radial velocity and computed the ensuing evolutions. 
Not only axisymmetric but non-axisymmetric perturbations have been also im- 
posed. We have applied mode analysis to the non-spherical deformation of the 
shock surface, using the spherical harmonics. We have found that (1) the growth 
rates of SASI are degenerate with respect to the azimuthal index m of the spher- 
ical harmonics Y™, just as expected for a spherically symmetric background, (2) 
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nonlinear mode couplings produce only m = modes for the axisymmetric per- 
turbations, whereas m 7^ modes are also generated in the non-axisymmetric 
cases according to the selection rule for the quadratic couplings, (3) the nonlinear 
saturation level of each mode is lower in general for 3D than for 2D because a 
larger number of modes are contributing to turbulence in 3D, (4) low I modes 
are dominant in the nonlinear phase, (5) the equi-partition is nearly established 
among different m modes in the nonlinear phase, (6) the spectra with respect to 
I obey power laws with a slope slightly steeper for 3D, and (7) although these 
features are common to the models with and without a shock revival at the end 
of simulation, the dominance of low I modes is more remarkable in the models 
with a shock revival. 

Subject headings: supernovae: collapse — neutrinos — hydrodynamics — insta- 
bility 



Introduction 



Many effo rts have been ma d e for t he multi-dime n sional modeling of core-collapse su- 
pernovae (see IWoosley fc Jankal (120051 ); iKotake et al.l (120061 ) for reviews), urged by accu- 
mulated observationajj^dejjces that c ore-collapse supernovae are globally aspherical com- 
monly (jWang et al.l Il996l . 1200 ll . |2002| ) . Vari ous mechanisms to p roduce the asymmetr y 



have been discussed so far: convection (e.g., iHerant et al.l (119941 ); 



Janka fc Mueller! (119961 ) ) . magnetic field and rapid rotation (see, e.g.. lKotake et al.l (120061 ) for 



Burrows et a 



4l995|); 



collective references) , standing (/stationary /spherical) accretion shock instability, or SASI , 



(IBlondin et al. 




2003 



Foglizzo et al. 



Scheck et all 120041 : IBlondin fe Mezzacappal I2OO6I : lohnishi et all bood. 



200 



6), and g-mode oscillations of protoneutron stars ( Burrows et al. 



Note, however, that most of them have been investigated with two-dimensional (2D) 



simulations. 



Recently a 3D study on SASI was reported by IBlondin &; Mezzacappal (120071 ). In 2D, 
the shock deformation by SASI is described with the Legendre polynomials Pi{0), or the 
spherical harmonics ^ m (#, <p) with m = 0. Various numerical simulations have demonstrated 
unequivocally that the I = 1 mode is dominant and a bipolar sloshing of the standing shock 
wave occurs with pu l sational strong expansions and contraction s along the symmetry axis 



(IBlondin et al.l 120031 ; IScheck et a 



2004 



Blondin & Mezzacappal 120061 : lOhnishi et al.l 12006 



Blondin fc Mezzacappa J2007h observed the dominance of 



20071 ). In 3D, on the other hand 
a non-axisymmetric mode with I — 1, m — 1, which produces a single- armed spiral in the 
later nonlinear phase. They claimed that this "spiral SASI" generates a rotational flow in 
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the accretion (see also iBlondin fc Shawl (120071 ) for 2D computations in the equatorial plane) 
and that it may be an origin of pulsar spins. 

There are many questions on 3D SASI remaining to be answered yet, however. For 
example, we are interested in how the growth of SASI differs between 3D and 2D among 
other things. In particular, the change in the saturation properties should be made clear. 
This will be the focus of th is paper. Another issue will be the generation of rotation in 
the accretion flow by SASI (IBlondin fc Mezzacappal 120071 ). Its efficiency and possible cor- 
relation with the net linear momentum should be studied more in detail and will be the 



subject of our forthcoming paper (llwakami et al.l 120071 ). Incidentally, it is noted that the 
neutrino heating and cooling we r e ent irely ignored and the flow was assumed to be isen- 
tropic in IBlondin fc Mezzacappal (120071 ). but the implementation of these physics is helpful 
for considering the impli cations not only for the shock revival but also for the nucleosynthe- 
sis fcifonidis et al.ll2006h . 



In this pape r, we have perfor med 3D hydrodynamic simulations, employing a realistic 
equation of state (jShen et al.lll998l ) and taking into account the heating and co oling of matter 



via n e utrino emissions and absorptions on nucleons as done in our 2D studies (lOhnishi et al. 



20061 . 120071 ) . The inclusion of the neutrino heating allows us to discuss how the critical 
luminosity for SASI-triggered explosion could be changed in 3D from those in 2D. To answer 
the questions we raised above, we have varied the initial perturbations as well as the neutrino 
luminosity and compared the growth of SASI between 2D and 3D in detail, conducting a 
mode analysis not only for the linear phase but also for the nonlinear saturation phase. 

The plan of this paper is as follows: In Section [2], we describe the models and numerical 
methods. We show the main numerical results in Section [3j We conclude the paper in 
Section |U 



2. Numerical Models 



T he numerical meth ods employed in the present paper are based on the code ZEUS- 
MP/2 (IHayes et al.l 120061 ) . which is a computational fluid dynamics code for the simulation 
of astrophysical phenomena, parallelized by the MPI (message-passing) library. The ZEUS- 
MP/2 code employs the Eulerian hydrodynamics algorithms based on the finite-difference 
method with a staggered mesh. In this study, we have modified t he original code substantial ly 
according to the prescriptions in our preceding 2D simulations (lOhnishi et al.ll2006l . 120071 ) . 



We consider spherical coordinates (r, 9, <fi) with the origin at the center of the proto- 
neutron star. The basic evolution equations describing accretion flows of matter attracted 
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by a proto-neutron star and irradiated by neutrinos emitted from the proto-neutron star are 
written as follows, 

dp 



P 



dv 

dt 



dt 

-- -VP 



+ p V • v = 0, 

p V $ - V ■ Q, 



d fe 

p dAp 



-PV ■ v + Q E - Q : Vv, 



dY P 



dt 



- Qn, 
GM in 



(1) 
(2) 
(3) 
(4) 
(5) 



where p, v, e, P, Y e , and $ are the density, velocity, internal energy, pressure, electron fraction, 
and gravitational potential, respectively. G is the gravitational constant. The self-gravity of 
matter in the accretion flow is ignored. Q is the artificial viscous tensor. Q E and Qn represent 
the heating/cooling and electron source/sink via neutrino absorptions and emissions by free 
nucleons. The Lagrangian derivative is denoted by d/dt = d/dt + v ■ V. The tabula ted 
realistic equation of state based on the relat ivistic mean field th eory (IShen et al.lll998l ) is 
implemented according to the prescription in lKotake et al.l (120031 ) . The mass accretion rate 
and the mass of the central object are fixed to be M = 1 M & s _1 and M- m = 1.4M , 
respectively. The neutrino heating is estimated under the assumptions that neutrinos are 
emitted isotropically from the cent ral object and that the neutrino flux is not affected by 
local absorptions and emissions (See Ohnishi et al. ( 20061 )). We consider only the interactions 
of electron-type neutrinos and anti-neutrinos. Their temperatures are also assumed to be 
constant and set to be T Uc = 4 MeV and T Pc = 5 MeV, which are the typical values in the post- 
bounce phase. The neutrino luminosity is varied in the range of L v = 6.0 — 6.8 x 10 52 ergs s -1 . 

The spherical polar coordinates are adopted. In the radial direction, the computational 
mesh is nonuniform, while the grid points are equally spaced in other directions. We use 
300 radial mesh points to cover r in < r < r out , where r in ~ 50 km is the radius of the inner 
boundary located roughly at the neutrino sphere and r out = 2000 km is the radius of the 
outer boundary, at which the flow is supersonic. 30 polar and 60 azimuthal mesh points are 
used to cover the whole solid angle. In order to see if this angular resolution is sufficient, 
we have computed a model with the 300 x 60 x 120 mesh points and compared it to the 
counterpart with the 300 x 30 x 60 mesh points. As shown later in Appendix [Bj the results 
agree reasonably well with each other both in the linear and nonlinear phases. Although the 
computational cost does not allow us to do the convergence test further, we think that the 
resolution of this study is good enough. 



- 5 - 



We use an artificial viscosity of tensor-type described in Appendix |X] instead of the 
von Neumann and Richtmyer-type that was originally employed in ZEUS-MP/2. For 3D 
simulations with the spherical polar mesh, we find it far better to employ the former than the 
latter to prevent the so-called carbuncle instability (IQuirk 1119941 ). which we observe around 
the shock front near the symmetry axis, 9 ~ 0, n. With the original artifitial viscosity, an 
appropriate dissipation is not obtained in the azimuthal direction fo r the shear flow resultin g 
from the converging accretion particularly when a fine mesh is used ( jStone fc Normanlll992l ). 
We have applied this artificial viscosity a lso to axisymmetric 2D simulations and reproduced 
the previous results (jOhnishi et al.ll2006l ). 



Figure [T] shows the radial distributions of various variables for the unperturbed flows. 
The spherically symmetr ic steady accre t ion fl ow through a standing shock wave is prepared 
in the same manner as in lOhnishi et al.l (120061 ). Behind the shock wave, the electron fraction 
is less than 0.5 owing to the electron capture, and a region of negative entropy gradient with 
positive net heating rates is formed for the neutrino luminosities, L v = 6.0— 6.8x 10 52 ergs s _1 . 
The values of these variables on the ghost mesh points at the outer boundary are fixed to 
be constant in time while on the ghost mesh points at the inner boundary they are set to be 
identical to those on the adjacent active mesh points, except for the radial velocity, which is 
fixed to the initial value both on the inner and outer boundaries. 

In order to induce the non-spherical instability, we have added a radial velocity per- 
turbation, 5v r (r,9,<p), to the steady spherically symmetric flow according to the following 
equation, 

v r (r,9,cj ) ) = v l r D (r) + 6v r (9, ( j ) ), (6) 

where vl D {r) is the unperturbed radial velocity. In this study, we consider three types of 
perturbations: (1) an axisymmetric (I — l,m — 0) single-mode perturbation, 



(7) 




Sv r (9,4>) oc \l^cos9-v 1 r D (r), 



(2) a non-axisymmetric perturbation with 1 = 1, 



6v r (r, 9, (f>) oc 




47T 



cos 6* 




8tt 



sin 9 cos « 



and (3) a random multimode perturbation, 

Sv r (9, cf)) oc rand ■ v]P{r) (0 < rand < 1), 



(8) 



(9) 



where "rand" is a pseudo-random number. The perturbation amplitude is set to be less than 
1% of the unperturbed velocity. It is noted that there is no distinction between m = 1 and 
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— 1 modes as long as the initial perturbation is added only to the radial velocity, as is the 
case in this paper. To put it another way, the m = ±1 modes contribute equally. Hence they 
are referred to as the \m\ — 1 mode in the following. On the other hand, differences do show 
up when the initial perturbation is added also to the non-radial vel ocity components, which 
case is important in discussing the origin of pulsar spins proposed bylBlondin fc Mezzacappa 



20071 ). All the models 



(120071 ) and will be the subject of our forthcoming paper (jlwakami et al 
presented in this paper are summarized in table HJ 

In the next section, we perform the mode analysis as follows. The deformation of the 
shock surface can be expanded as a linear combination of the spherical harmonics components 
17"(M): 

oo I 

fls(M) = EE WiW), (io) 

1=0 m=-l 

where Y™ is expressed by the associated Legendre polynomial P™ and a constant Kf 1 given 
as 

(11) 



yra = #™pm( cos 0) gim* 



'2/ + 1 (l-m)\ 



4vr (l + m)l 

The expansion coefficients can be obtained as follows, 



(12) 



m 
C l 



2lT /*7T 



o jo 



d9 sin R s {6, <f>) Yr* (0, <p), 



(13) 



where the superscript * denotes complex conjugation. 



3. Results and Discussions 

3.1. Axisymmetric Single-Mode Perturbation 

Before presenting the results of 3D simulations, we first demonstrate the validity of our 
newly developed 3D code. We compare the axisymmetric flows obtained by 2D and 3D 
simulations for the axisymmetric (/ = l,m = 0) single-mode perturbation. By 2D simu- 
lations we mean that axisymmetry is assumed and computations are done in the meridian 
section with all 0- derivatives being dropped, whereas in 3D simulations we retain all these 
derivatives and 3D mesh is employed. This validation is important because numerical errors 
may induce azimuthal motions even for the axisymmetric initial conditions. Hence we have 
to confirm that azimuthal errors do not appear in the non-axisymmetric simulation. 
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Figure [2] shows the time evolutions of the average shock radius Rs- The 2D and 3D 
results are displayed in Figure [2] (a) and (b), respectively. The average shock radius Rs 
is obtained from the expansion coefficient c° in Eq. (fTUj) multiplied by K®. The solid, 
broken, and dotted lines correspond to the neutrino luminosities L v of 6.0, 6.4 and 6.8 x 
10 52 ergs s _1 , respectively. It should be first mentioned that one cannot expect an perfect 
agreement between two computations of the exponential growth of the instability followed 
by the turbulent motions through mode couplings as considered here. It is still obvious that 
the results of the 2D and 3D simulations agree on essential features. In particular, both in 
the 2D and 3D results, we see that the Rs is settled to an almost constant value for L v = 6.0 
and 6.4 x 10 52 ergs s _1 , whereas it continues to increase for L v = 6.8 x 10 52 ergs s^ 1 . 

Figure [3] presents the normalized amplitudes |cJ™/co| as a function of time for Model I 
(L„ = 6.0 x 10 52 ergs s" 1 ). The 2D and 3D results are displayed in Figure [3] (a) and 
(b), respectively. The red, blue, black, and gray solid lines correspond to the modes of 
(l,m) = (1,0), (2,0), (3,0) and (4,0), respectively. As can be seen, the time evolution can 
be divided into two phases. One is the linear growth phase, in which the amplitude of the 
initially imposed mode with (/, m) = (1, 0) grows exponentially. It lasts for ~ 100 ms. It 
is also found that higher modes are generated by mode couplings and grow exponentially 
during this phase. Then starts the second phase, which is characterized by the saturation 
of amplitudes of the order of 0.1. In this phase, the accretion flow becomes turbulent. It is 
interesting to note that in this nonlinear saturation phase the amplitude of the mode with 
(I, m) = (2, 0) is almost of the same order as the initially impo sed mode with (/, m) = (1, 0) 



and is dominant over other modes, which fact was observed in lOhnishi et al.l (120061 ) 



Most important of all, however, is the fact that all the m/fl modes are not produced, 
implying that the perturbed flow retains axisymmetry, which is a necessary condition for 
the numerical study of non-axisymmetric instability. Although the results of the 2D and 3D 
simulations are not identical again, the essential features such as the linear growth rate of the 
initial perturbation (I = l,m = 0) and the production of other modes via nonlinear mode 
couplings as well as the saturation levels in the nonlinear phase are all in good agreement for 
the two cases. The quantitative differences between the 2D and 3D results mainly originate 
from the difference in the time steps, which depend on the minimum grid width. Note that 
in addition to the 300 radial and 30 polar mesh points used for the 2D computations, 60 
azimuthal mesh points are employed in the 3D simulatioins and, as a result, smaller time 
steps are usually taken for the 3D computations. 
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3.2. Non-Axisymmetric Perturbation with / = 1 

Now we discuss the results of genuinely 3D simulations, in which the non-axisymmetric 
perturbation with (I = 1, \m\ = 1) is added to the axisymmetric (I — 1, m = 0) perturbation. 

Figure H] shows the time evolutions of some of the normalized amplitudes |c™/cq| for 
Model II with the neutrino luminosity L u = 6.0 x 10 52 ergs s -1 . The red, yellow, blue, light 
blue, green, black, brown, violet, and pink solid lines denote the amplitudes of the modes with 
(I, \m\) = (1, 0), (1, 1), (2, 0), (2, 1), (2, 2), (3, 0), (3, 1), (3, 2) and (3, 3), respectively. In the 
linear phase, the initially imposed modes with (/, |m|) = (1, 0) and (1, 1) grow exponentially 
just as in the axisymmetric single-mode perturbation. It should be noted that the growth 
rate of (I, \m\) = (1,1) mode is identical to that of (l,m) = (1,0) mode. This is just as 
expected for the spherically symmetric background, for which the modes with the same I 
but different m are degenerate in the linear eigen frequency. 

It is observed in the figure that the modes with (I, \m\) = (2,0), (2, 1), (2,2) are first 
produced by the nonlinear mode-couplings of (I, |m|) = (1,0), (1, 1). Then the modes with 
(/, \m\) = (3, 0), (3, 1), (3, 2), (3, 3) can be produced by the couplings of the first-order modes 
with (l,\m\) = (1,0), (1,1) and the second-order modes with (l,\m\) = (2, 0), (2, 1), (2, 2). 
Even higher order modes are produced subsequently toward the non-linear saturation but 
are not shown here. Although the branching ratios should be investigated more in detail 
before giving any conclusion, the above sequence of the mode generations strongly suggests 
that the nonlinear coupling is mainly of quadratic nature. 

In the nonlinear phase that begins at t ~ 200ms, these mode amplitudes are saturated 
and settled into a quasi-steady state. As in the axisymmetric case, the I = 1 modes both with 
m = and |m| = 1 are dominant over other modes in this stage for the non-axisymmetric 
case. One thing to be noted here, however, is the fact that the saturated amplitudes for 
Model II are lower than those for Model I in general (see Figure [3] (b)). This is, we think, 
because the turbulent energy is shared by a larger number of modes in the non-axisymmetric 
case than in the axistmmetric case. To demonstrate this more clearly, we have added the 
random perturbation to the radial velocity in the axisymmetric Model I at t — 400 ms, by 
which time the axisymmetric nonlinear turbulence has fully grown. We refer to this model 
as Model HI. The entire time evolutions of the normalized amplitudes |c™/c(j| for Model HI 
are shown in Figure It is clear that the axisymmetric m = mode amplitudes are reduced 
after the additional perturbation is given and the non-axisymmetric m^0 modes grow to 
the saturation level at t ~ 450 ms. The power spectra of the turbulent motions will be 
discussed more in detail later. 
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3.3. Random Multi-Mode Perturbations 

3.3.1. Dynamical Features and Critical Luminosity 

Having understood the elementary processes of the linear growth and nonlinear mode 
couplings in the previous section, we will now discuss the 3D SASI induced by the random 
multimode perturbations, which are supposed to be closer to the reality. In this subsection, 
we will show the typical dynamics, paying attention to the time evolution of the shock wave 
and, hence, to the critical luminosity, at which the stalled shock is revived. 

Figure [6] shows the time variations of the average shock radius Rs for Models IV, V and 
VI with the neutrino luminosities, L v = 6.0,6.4 and 6.8 x 10 52 ergs s _1 , respectively. It is 
found that for L v = 6.0 and 6.4 x 10 52 ergs s -1 , the shock is settled to a quasi-steady state 
after the linear growth whereas it continues to expand for L u = 6.8 x 10 52 ergs s _1 , which 
is also true for the axisymmetric counterpart displayed in the right panel of the figure for 
comparison. This implies that the critical neutrino luminosity is not much affected by the 
change from 2D to 3D. In the following we refer to the models with L v = 6.0 and 6.4 x 10 52 
ergs s _1 as the non-explosion models and the model with L v = 6.8 x 10 52 ergs s _1 as the 
explosion model and look into their dynamical features in turn. The mode analyses will be 
done in sections 13.3.21 and 13.3.31 

Figure [7] shows the sideviews of the iso-entropy surfaces together with the velocity 
vectors in the meridian section at four different times for the non-explosion model with 
L v = 6.0 x 10 52 ergs s -1 (Model IV). The hemispheres (0 < < 7r) of eight iso-entropy 
surfaces are superimposed on one another and the outermost surface almost corresponds 
to the shock front. The initial perturbations grow exponentially in the linear phase ((a) 
t = 40 ms). At the end of the linear phase, a blob of high entropy is formed and grows 
subsequently ((b) t = 93 ms). High entropy blobs are continuously generated and the 
nonlinear phase begins ((c) t = 100 ms). Circulating flows occur inside the blobs while 
high-velocity accretions onto a proto-neutron star surround the blobs. These blobs expand 
and shrink repeatedly, being distorted, splitted and merged with each other inside the shock 
wave ((d) t = 350 ms). Reflecting these complex motions, the shock surface oscillates in all 
directions but remains alomost spherical for the non-axisymmetric model. This is in sharp 
contrast to the axisymmetric ca se, in which large-am p litude oscillations occ ur mainly in the 



direction of the symmetry axis (IBlondin et al.l (120031 ); lOhnishi et al.l (120061 )) 



Figure [8] displays the snapshots of the iso-entropy surfaces and the velocity vectors for 
the explosion model with L u = 6.8 x 10 52 ergs s -1 (Model VI). As in the non-explosion 
model, the sequence of events starts with the linear growth of the initial perturbation inside 
the shock wave ((a) t = 40 ms). In the explosion model, however, many high entropy blobs 
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emerge simultaneously much earlier on ((b) t = 70 ms) than in the non-explosion model. 
Then, these blobs repeat expansions and contractions, being distorted, splitted and merged 
together just as in the non-explosion counterpart ((c) t — 80 ms). After that, some of the 
blobs get bigger, absorbing other blobs ((d) t = 350 ms) and, as a result, the shock radius 
increases almost continuously up to the end of the computation (t = 400 ms) as already 
demonstrated in Figure [6] (a). At this point of time, the shock surface looks ellipsoidal 
rather than spherical. However, the major axis is not necessarily aligned with the symmetry 
axis and the flow is not symmetric with respect to this major axis, either. 

3.3.2. Mode Spectra 

Next we look into the spectral intensity more in detail. As a standard case, we take 
Model IV with L v = 6.0 x 10 52 ergs s _1 . Figure [9] shows the time evolutions of the normalized 
amplitudes |c™/cq| with m = and compares them with the axisymmetric counterparts in 
Model VII. Note that the m ^ modes also exist in the non-axisymmetric model, which 
are not shown in the figure but will be discussed in the following paragraphs. As can be 
clearly seen, the amplitude of each mode grows exponentially until ~ 100 ms, which is the 
linear phase. Note in particular that the growth rate and oscillation frequency for the I = 1 
mode are the same as those obtained for the model with the single mode perturbation in 
the previous section. After ~ 100 ms, the mode amplitudes are saturated and the evolution 
enters the nonlinear phase with a clear dominance of the modes with I — 1, 2. It is also 
evident in the figure that the saturation level is lower in general for the non-axisymmetric 
case than for the axisymmetric one, which confirms the claim in the previous section based 
on the results for the single-mode perturbation. 

In Figure [T0l we display the snapshots at four different times of the spectral distributions 
both for the non-axisymmetric (the left panels) and axisymmetric (the right panels) cases. 
The uppermost panels correspond to the linear phase and the intensity is distributed rather 
uniformly over all modes. As the time passes, however, the amplitudes of low / modes grow 
much more rapidly than those with higher /'s (the second and third rows in the figure). One 
can see a similarity in the time evolutions between the non-axisymmetric and axisymmetric 
models. After the nonlinear phase starts (the lowermost panels), the growths of all modes 
are saturated and the spectral distributions are settled to be quasi-steady. It is again obvious 
in these panels that the low / modes are dominant in the nonlinear phase and the saturation 
level is lower in the non-axisymmetric case. 

Now we turn our attention to the quasi-steady turbulence in the nonlinear phase. Shown 
as a function of I in Figure [IT] are the power spectra |c™/c(]| 2 averaged over the interval from 
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t = 150 ms to 400 ms. Both the non-axisymmetric and axisymmetric cases are displayed by 
the symbols, o and •, respectively. In the left panel, modes with different m's are plotted 
separately whereas they are summed up in the right panel. It is found that the time-averaged 
power spectra are not much different among the modes with the same I but with a different 
m. This implies that the equi-partition of the turbulence energy is nearly established among 
the modes with the same I. This will have an import ant ramification for t he origin of pulsar 



spin and will be discussed in our forthcoming paper (jlwakami et al.l 120071 ). 



The right panel of the figure demonstrates that the time-averaged power spectrum 
summed over m for the non-axisymmetric case is not much different from that for the ax- 
isymmetric counterpart. This means that the total turbulence energy does not differ very 
much between the axisymmetric and non-axisymmetric cases. As a result, the power of each 
mode in the non-axisymmetric case is smaller roughly by a factor of 21 + 1 than that of the 
same I mode in the axsymmetric case. This is the reason why we have observed that the 
saturation level of the nonlinear SASI is smaller in the non-axisymmetric case than in the 
axisymmetric case and the shock surface oscillates in all directions with smaller amplitudes 
in the non-axisymmetric flow whereas it sloshes in the direction of the symmetry axis with 
larger amplitudes in the axisymmetric case. The difference in the fluctuations of the average 
shock radius seen in Figure [6] (a) and (b) is also explained in the same manner. 

Another interesting feature observed in the figure is the fact that the time-averaged 
power spectra obey a power-law at I > 10 for both the non-axisymmetric and axisymmetric 
models. Two straight lines in Figure fTTl (a) are the fits to the data for I > 10, each obtained 
for the axisimmetric and non-axisymmetric models. The powers are found to be —5.7 and 
—4.3 for the non-axisymmetric and axisymmetric cases, respectively. Although the difference 
is almost unity, which originates from the difference in multiplicity of the modes with the 
same I, the slope excluding this effect is still a little bit steeper in the non-axisymmetric case 
as seen in Figure [TT] (b). At the moment we do not know if this difference in slope is real or 
not and the power itself is also remaining to be explained theoretically. 



3.3.3. The Explosion Models 

So far we have been discussing the non-explosion models, in which the SASI is saturated 
in the non-linear phase and is settled to a quasi-steady state. In considering the possible 
consequences of SASI after the supernova explosion occurs, however, it is also interesting to 
investigate the explosion models, in which a shock revival occurs as a result of SASI. 

Figure [T21 shows the time evolutions of the normalized amplitudes |c™/cq| for the ex- 



- 12 - 



plosion model (Model VI) and compares it with the axisymmetric counterpart (Model VIII). 
The feature of the explosion models common to the axisymmetric and non-axisymmetric 
cases is as follows: The nonlinear stage is devided into two phases. The earlier phase is quite 
similar to the nonlinear phase that we have seen so far for the non-explosion models. The 
later phase, on the other hand, has a distinctive feature that the I = 1 mode becomes much 
more remarkable than other modes and the oscillation period tends to be longer as the shock 
radius gets larger and the advection-acoustic cycle, supposedly an underlying mechanism of 
the instability, takes longer. The later prominence of / = 1 mode is intriguing and is impor- 
tant to discuss the kick velocity of pulsar quantitatively. The theoretical account, however, 
is yet to be given. The differences between the non-axisymmetric and axisymmetric models, 
such as the saturation level, are just taken over to the explosion models. 

The power spectra averaged over the time intervals of 100 < t < 200 ms and 250 < 
t < 400 ms are presented in Figure [13] (a) and (b), respectively. As in the non-explosion 
models, the equi-partition among different m modes is again established for the explosion 
models. This is true even in the later nonlinear phase as seen in the right panel of the 
figure. As a result, the saturation level is smaller for the non-axisymmetric case than for 
the axisymmetric case as mentioned just above. The power spectra in the earlier nonlinear 
phase look much similar to those found in the non-explosion models, with the power law 
being satisfied for / > 10. In the later nonlinear phase, on the other hand, the power law is 
extended to much lower Z's. This is related to the late-time prominence of the / = 1 mode 
as mentioned above. In fact, the amplitudes in the lower I portion of the power spectrum 
are enhanced in 250 < t < 400. The mechanism of this amplification is remaining to be 
understood but might be related with the volume-filling thermal co nvection advocated fo r 



the late stage of the convective instability in the supernova core by iKifonidis et al.l (120061 ). 
As a result of this evolution of the spectrum, the shock tends to be ellipsoidal as the shock 
expands. 



3.4. Neutrino Heating 

The SASI is supposed to be an important ingredient not only for the pulsar's proper 
motions and spins but also for the explosion mechanism itself. In this section, we look into 
the neutrino heating in the non-axisymmetric SASI. 

We show in Figure [TH the color contours of the net heating rate in the meridian section 
for the non-explosion model with L v = 6.0 x 10 52 ergs s _1 (non-axisymmetric Model IV 
and axisymmetric Model VH) as well as for the non-axisymmetric explosion model with 
L u = 6.8 x 10 52 ergs s _1 (Model VI ). In the early phase, the cooling region with negative 
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net heating rates exists around the proto neutron star while the heating region extends over 
the cooling region up to the stalled shock wave in all the cases. As the time passes and the 
SASI grows, a pocket of regions with positive but relatively low net heating rates appear. 
These regions correspond to the high entropy blobs (high entropy rings for the axisymmetric 
case), where the circulating flow exists observed in Figure [7J (b), (c), and (d). Since the 
neutrino emission in these regions is more efficient than the surroundings, the net heating 
rate is rather low. 

Although the critical neutrino luminosties are not much different between the non- 
axisymmetric and axisymmetric cases, the spatial distributions of neutrino heating are dif- 
ferent. In the axisymmetric case the shock wave oscillates up and down whereas it moves 
in all directions in the non-axisymmetric case. The oscillation amplitudes are larger in the 
former than in the latter in general as repeated mentioned. Reflecting this difference in 
the shock motions, the neutrino heating is enhanced in mainly in the polar regions in the 
axisymmetric case while in the non-axisymmetric case the heating rate is affected by SASI 
chiefly around the high entropy blobs. These are clearly seen in Figure [TH (see also Figure [7]). 

As described in section 13.3. 1\ the generation of the high entropy blobs is started around 
the end of the linear phase and continues during the nonlinear phase. Although the blobs 
repeat expansions and contractions, being sometimes merged or splitted, the turbulent mo- 
tions together with the neutrino heating are finally settled to a quasi-steady state in the 
non-explosion model. For the explosion model, on the other hand, the number and volume 
of high entropy blobs are increased much more quickly and, as a result, the heating region 
prevails, pushing the shock wave outwards and narrowing the cooling region near the proto 
neutron star. 

Figure [TBI shows the time evolutions of the net heating rates integrated over the region 
inside the shock wave. For the non-explosion model (Model IV, solid line), the volume- 
integrated heating rate grows until t ~ 150ms but is saturated thereafter, whereas it increases 
continuously for the explosion model (Model VI, dashed-dotted line). These behaviors of the 
volume-integrated heating rate are in accord with the time evolutions of the average shock 
radius as shown in Figure [6] (a)- It is clear from the figure that the heating rates are mainly 
affected by SASI during the nonlinear phase. For comparison, the spherically symmetric 
counterparts, Models DC with L v = 6.0 x 10 52 ergs s _1 and X with L v = 6.8 x 10 52 ergs s _1 
are also presented as dashed and dotted lines, respectively. Note that the last model does 
not lead to the shock revival even with this high neutrino luminosity. It is found that the 
SASI enhances the neutrino heating both for the explosion and non-explosion models. 
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4. Conclusion 

In this paper we have studied the non-axisymmetric SASI by 3D hydrodynamical sim- 
ulations, taking into account the realistic EOS and neutrino-heating and cooling. We have 
added various non-axisymmetric perturbations to spherically symmetric steady flows that ac- 
crete through a standing shock wave onto a proto-neutron star, being irradiated by neutrinos 
emitted from the proto-neutron star. The mode analysis has been done for the deformation 
of the shock surface by the spherical harmonics expansion. After confirming that our 3D 
code is able to reproduc e for the axisymmetr ic perturbations the previous results on 2D 



SASI that we obtained in lOhnishi et al.l (120061 ). we have done genuinely 3D simulations and 
found the followings. 

First, the model of the initial perturbation with (I, \m\) = (1, 0) and (1, 1) have demon- 
strated that the non-axisymmetric SASI proceeds much in the same manner as the axisym- 
metric SASI: the linear phase, in which the initial perturbation grows exponentially, lasts for 
about 100ms and is followed by the nonlinear phase, in which various modes are produced 
by nonlinear mode couplings and their amplitudes are saturated. It has been found that 
the critical neutrino luminosity, above which the shock revival occurs, is not much different 
between 2D and 3D. For the neutrino luminosities lower than the critical value, the SASI is 
settled to a quasi-steady turbulence. We have found that the saturation level of each mode 
in the non-axisymmetric SASI is lower in general than that of the axisymmetric counterpart. 
This is mainly due to the fact that the number of the modes contributing to the turbulence 
is larger for the non-axisymmetric case. The sequence of the mode generation, on the other 
hand, strongly suggests that the nonlinear mode coupling is chiefly of quadratic nature. 

Second, the simulations with the random multi-mode perturbations being imposed ini- 
tially have shown that the dynamics in the linear phase is essentially a superimposition of 
those of single-modes. Toward the end of the linear phase, high entropy blobs are generated 
continuously and grow, starting the nonlinear phase. We have observed that these blobs 
repeat expansions and contrations, being merged and splitted from time to time. In the 
non-explosion models, these nonlinear dynamics lead to the saturation of mode amplitudes 
and the quasi-steady turbulence. For the explosion models, on the other hand, the pro- 
duction of blobs proceeds much more quickly, followed by the oligarchic evolution, with a 
relatively small number of great blobs absorbing smaller ones, and as a result the shock 
radius increases almost monotonically. The spectral analysis has clearly demonstrated that 
low I modes are predominant in the nonlinear phase just as in the axisymmetric case. We 
have also shown that the equi-partition is nearly established among different m modes in 
the quasi-steady turbulence and that the spectrum summed over m in the non-axisymmetric 
case is quite similar to the axisymmetric counterpart. This implies that the larger number 
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of modes in the non-axisymmetric case is the main reason why the amplitude of each mode 
is smaller in 3D than in 2D. The power spectrum is approximated by the power law for 
I > 10. Although the slope is slightly steeper for the non-axisymmetric models, whether the 
difference is significant or not is unknown at present. 

We have seen in the explosion models, on the other hand, that the oscillation period of 
each mode becomes longer in the late nonlinear phase as the shock radius gets larger and the 
advection-acoustic cycle becomes longer. What is more interesting is the fact that in this 
late phase the dominance of low / modes becomes even more remarkable. Although this may 
be related with the volume-filling thermal convection, the mechanism is yet to be revealed. 
This spectral evolution leads to the global deformation of the shock surface to an ellipsoidal 
shape, whose major axis is not necessarily aligned with the symmetry axis, however. 

Finally, we have presented the neutrino heating in the 3D SASI. It has been seen that 
the volume-integrated heating rate is affected mainly in the nonlinear phase. The comparison 
with the spherically symmetric counterparts has confirmed that the SASI enhances the neu- 
trino heating also in the non-axisymmetric case. Although the critical neutrino luminosity 
in the non-axisymmetric SASI is not much changed from that for the axisymmetric case, the 
spatial distribution of neutrino heating is different in the non-linear phase. Relatively nar- 
row regions surrounding high entropy blobs are efficiently heated for the non-axisymmetric 
case, while wider regions near the symmetry axis are heated strongly in accord with the 
sloshing motion of shock wave along the symmetry axis for the axisymmetric case. For the 
non-explosion models, the high entropy blobs produced by neutrino heatings occupy the 
inside of shock wave, repeating expansion and contraction and being splitted and merged 
intermittently, but the flow is finally settled to a quasi-static state. For the explosion mod- 
els, on the other hand, the high entropy blobs are generated much more quickly and extend 
the heating region, pushing the shock outwards and narrowing the cooling region near the 
proto-neutron star. 

In the present study we have not obse rved a persistent segreg a tion o f angular momen- 



tum in the accretion flow such as found by lBlondin &: Mezzacappal (120071 ); iBlondin fc Shaw 



(120071 )) although instantaneous spiral flows are frequently seen. As discussed above, the 
equi-partition is nearly established among different m modes in our models. It should be em- 
phasized here again, however, that we have added the initial perturbations only to the radial 
velocity in this study and, as a result, the modes with m = ±1 are equally contributing. We 
will defer the analysis of the models with the initial non-axisymmetric pertu rbations being 



adde d also to the azimuthal component of velocity to our forthcoming paper (jlwakami et al. 



20071 ). in which we will also discuss a possible correlation between the kick velocity and spin 



of neutron stars if they are produced by the 3D SASI indeed. A lot of issues on SASI are 
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still remaining to be studied. In particular, the influences of rotation and magnetic field are 
among the top priorities and will be addressed elsewhere in the near future. 
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A. Tensor Artificial Viscosity 



We present the tensor-type artificial viscosity, Q, that we use in this paper. It is a direct ex- 
tention to 3D of the one employed bvlStone &: Norman! (119921 ) for 2D simulations. Following 
the notations in lStone &: Norman! (119921 ). Q is written as follows: 



Q 



'l 2 pW ■ v [V« - | (V • v) e] if V ■ v < 







otherwise 



(Al) 



where / denotes a constant with a dimension of length, Vv = (vij+Vj A )/2 is the symmetrized 
velocity-gradient tensor, e is the unit tensor. Dropping the off-diagonal components, we 
utilize only the diagonal components of Q in this study, which are written in the finite 
difference form as 



Q 11 i,j,k = iUkdi,j,k (v • v) ijk \iy V{ll) ). hk - |(v • v\ jth 

Q22 itj>k = ll jjk d i>j>k (V • v) lhk [(V« (a2) ) . - | (V • v\. tk 
Q33 UfJfc = ll hk d hhk (V • v) i)j>k \(Vv m ) - I (V • v) i; . k 



(A2) 
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where dij t k stands for the density at the site specified by (i,j, k) and Uj )k is written as 

lijjt = max(C2 dxlai,C% dxla^^C^ dx3a k ). (A3) 

Here C 2 is a dimensionless constant controlling the number of grid points, over which a shock 
is spread, and dxla iy dx2aj, and dx3a k are the grid widths of the "a-mesh" at the i-, j-, and 
fc-th grid points in the r, 9, and <fi directions, respectively. The "a-mesh" and "b-mesh" are 
disti nguished in ZEUS-MP/ 2 and are defined on the cell-edge and cell-center, respectively 



see 



Stone fc Norman! (119921 ) for more details). 



The artificial viscous stress —V - Q and artificial dissipation — Q : Vv in the momentum 
equation (j2J and energy equation (j3j) are given as follows, 



d_ 

1 9 I 2 r\ \ 1 d 932 



(V • Q) (2) = -^rj- + (A4) 

1 ' #2#32 0x2 929m 0x2 

(v-Q) (3) = — fa, 

1 ' 92932 dx 3 



Q : X7v = l 2 P V ■ v\ (Vv(n) - Vv (22) ) 2 + (VV(u) - Vv (33 )) 2 + (Vv (33 ) - Vv (22 )) 2 , (A5) 



where g 2 = r, g%\ = r, and = sin 9 a re the metric componen ts for the spherical co- 
ordinates, (xi,X2,Xs) = (r,9,<fi) (see again 
discretized as 



Stone fc Norman! (11992I )). These equations are 



g2blg31biQll i>jjk - g2bl l9 31b i . l Qll 



i—l,j,k 



g2ajg31aidxlbi 



ry Q) g32b)Q2\ h k - g3jg zl g22 M _i i fc Ql^+Ql^,^ f dg32aA 

[ ' ^ 2 «> g 2b l g32a 2 dx2b j 2g2b i g32a j \ dx 2 ) 



(X7 c\\ - Q33ij.k Q33ij : k—i 
I • ^)3,i,j,k ~ g 2b il g 32b j dx3b k ' 



(Q : V«) M - fc =l 2 i jik dij, k {V ■ «)ij,*|[[(V«( 1 i))ij, fc - (Vv (22) )i Jjfe ] 2 

+ [(VV(n))ij,fc - (V«(33)),,i,fc] 2 + [(V«(33))ij,Jfc - (VV(22))i,i,fc] 2 ], 

where p2a, g31a, and g32a are <?3i> and #32 defined on the "a-mesh", whereas g2b, g31b, 
and g32b are those defined on the "b-mesh". dxlbi, dx2bj, and dx3b k represent the width of 
the "b-mesh" at the i-, j-, and fc-th grid points in the r, 9, and <fi directions, respectively. 
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Finally the velocity-gradient tensor (Vw(n))ij,it! (VV( 2 2))ij,fe, and (Vf^))^^ are given by 
the following equations: 



Ul 



* J,fc efccla. 
_' u 2jj +lifc — 



ifti,j,k , «lj j,fc + «li+ij,fc / dg2bj \ 



+ 



V 9xi / 



g2bidx2aj 2g2bi 
v\j,k+x ~ v3i,j, k . v2i j- k + v2 i>j+1>k f dg32bj\ vl i)jtk + vl i+1)jtk f dg31bA 



^ WJi&k - g 2b i g32b j dx3a k 2g31b i g32b j \ dx 2 



+ 



2g31h 



V 9xi ) 
{Ah 



B. Numerical Convergence Tests 



In order to see if the numerical resolution employed in the main body is sufficient, we 
increase the number of angular grid points and compare the results. Figure [16] (a) shows the 
time evolutions of the normalized amplitudes | cj" / Cq | in the linear phase. In this comparison, 
we impose the I = 1, \m\ — 1 perturbation initially. We refer to the model with 300 x 30 x 60 
mess points as MESHO and to that with 300 x 60 x 120 grid points as MESHlin the figure. 
It is found that the linear growth rates agree with each other reasonably well although 
the coarser mesh slightly overestimates the growth time. Figure [16] (b) presents the power 
spectra |c™/cq| 2 that are time-averaged over the nonlinar phase. The random perturbation 
is imposed in this case. It is again clear that the results for MESHO is in good agreement 
with those for MESH1. 

It should be mentioned that for MESHO it takes 32 parallel processors about 1.5 days to 
compute the evolution up to t = 400 ms, while MESH1 needs 22 days even for 128 parallel 
processors. This is partly because of the difference in the Courant numbers, which are set to 
be 0.5 for MESHO but to be 0.1 for MESH1 to better use the tensor- type artificial viscosity. 
Although this severe limit of CPU time does not allow us to do more thorough convergence 
tests, we think, based on the results of the tests shown above, that our results given in this 
paper are credible. 
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Table 1: Summary of all the models. 



Model 


Perturbation 


Nuetrino Luminosity L v 






[10 52 ergs s _1 ] 


I 


single- mode, I — l,m — 


6.0 


n 


multi-mode, I — 1, m — and I = 1, m = 1 


6.0 


m 


single- mode, I — l,m — with 


6.0 




random perturbation at t = 400 ms 




IV 


random perturbation 


6.0 


V 


random perturbation 


6.4 


VI 


random perturbation 


6.8 


vn 


axisymmetric random perturbation 


6.0 


vm 


axisymmetric random perturbation 


6.8 


DC 


none 


6.0 


X 


none 


6.8 
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Fig. 1. — Profiles of various variables for the unperturbed spherically symmetric accretion 
flows with different neutrino luminosities 
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(b) 3D axisymmetric perturbation 



Fig. 2. — The time evolutions of the average shock radius Rs in 2D (left panel) and 3D 
(right panel) for the axisymmetric I — l,m — single-mode perturbation. 
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(a) 2D axisymmetric perturbation (b) 3D axisymmetric perturbation 



Fig. 3. — The time evolutions of the normalized amplitudes |c™/cq| for Model I with the 
axisymmetric I — 1, m — single-mode perturbation in 2D (left panel) and 3D (right panel). 
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Fig. 4. — The time variations of the normalized amplitudes |cJ™/cq| for Model II, in which 
the modes with (/, |m|) = (1, 0), (1, 1) are initially imposed. 
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Fig. 5. — The time evolutions of the normalized amplitudes |c™/cq| for Model m, in which 
we have added the random perturbation to the radial velocity at t — 400ms. 
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(a) 3D non-axisymmetric random perturbation(b) 2D axisymmetric random perturbation 

Fig. 6. — The time evolutions of the average shock radius R$ for Models IV, V and VI, 
in which the random multi-mode perturbation is imposed in 3D (left panel) and 2D (right 
panel). 
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(c) t = 100 ms (d) t = 390 ms 

Fig. 7. — The iso-entropy surfaces and velocity vectors in the meridian section for Model IV. 
The hemispheres (0 < <fi < tt) of eight iso-entropy surfaces are superimposed on one another 
with the outermost surface nearly corresponding to the shock front. 
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(c) t = 80 ms (d) t = 350 ms 



Fig. 8. — The iso-entropy surfaces and velocity vectors in the meridian section for the 
explosion model (Model VI). Note that the displayed region is 2.5 times larger for panel (d). 
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(a) Model IV (b) Model MR 



Fig. 9. — The time evolutions of the normalized amplitudes |c™/cq| for Model IV (non- 
axisymmetric) and VH (axisymmetric). Note that the m ^ modes also exist in the non- 
axisymmetric model but are not shown in this figure. 
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Fig. 10. — The normalized amplitudes |c™/cq| for Model IV (left panels) and Model VII (right 
panels) at different times. Note that the time-averaged values are plotted in panels (d) and 
(h). 
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Fig. 11. — The time- averaged power spectra for Model IV and VII. The average is taken over 
t = 150 — 400 ms. In the left panel, modes with different m's are plotted separately whereas 
they are summed up in the right panel. 
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(a) Model VI (b) Model Vffl 

Fig. 12. — The time evolutions of the normalized amplitudes |c™/c[}| for the explosion models 
(Model VI (non-axisymmetric) in the left panel and Model VTH (axisymmetric) in the right 
panel. 
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(b) t = 250 - 400 ms 
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Fig. 13. — The time-averaged power spectra for the explosion models, Model VI (non- 
axisymmetric) and Vffl (axisymmetric). The averages are taken (a) for t = 100 — 200 ms 
and (b) for t = 250 — 400 ms, respectively. 
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(b) Model IV, t = 390 ms 
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(f) Model VI, t = 350 ms 



Fig. 14. — The color contours of the net heating rate in the meridian section for Model IV 
(top panels), Model VH (middle panels) and Model VI (bottom panels) at different times. 
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Fig. 15. — The time variations of the net heating rate integrated over the region inside the 
shock wave. Models IV and VI are compared with the spherically symmetric counterparts, 
Models DC and X. 
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(b) multimode random perturbation 



Fig. 16. — The convergence tests, (a) The time evolutions of the normalized amplitudes 
|c™/c[}| in the linear phase. The single-mode perturbation with I = 1, \m\ = 1 is imposed 
initially. The models with 300 x 30 x 60 and 300 x 60 x 120 mess points are referred to 
as MESHO and MESH1, respectively, (b) The time-averaged power spectra |c™/cq| 2 . The 
average is taken over 150 < t < 400 ms. In this comparison, the random multi-mode 
perturbation is imposed. 



